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SUMMARY 


Two new inverse boundary -layer computational procedures, which permit nonsingu- 
lar solutions at separation and reattachment, are presented. In the first technique, which 
is for incompressible flow, the displacement thickness is prescribed and the pressure is 
deduced implicitly along with the solution. Comparisons are made by using this pro- 
cedure with the turbulent experimental data of Chu and Young, and excellent agreement 
is obtained. The purpose of the second technique is to obtain an inverse computational 
procedure for compressible separated flows that is nonsingular at separation. In this 
formulation, a perturbation mass flow is prescribed for the compressible boundary- layer 
equations and the pressure is deduced implicitly along with the solution. Laminar and 
turbulent computations made by using this technique are presented and are typical of sep- 
arated flow. In both inverse procedures, finite -difference techniques are used along with 
Newton iteration. The present boundary-layer procedures can be used either in the 
inverse or direct (pressure prescribed) mode, depending on the boundary conditions; in 
either mode the solution procedure is no more complicated than conventional boundary- 
layer computations. Both inverse techniques presented herein appear to be well suited 
for complete viscous-inviscid interaction computations. 

INTRODUCTION 

In recent years, significant progress has been made in the solution of flow fields 
which contain regions of separated flow. Despite the diversity of approaches, which have 
been reported, there seems to be a clear distinction that most techniques are either based 
on the Navier-Stokes equations with relatively large computer costs or more approximate 
procedures in which the simultaneous solution of the boundary -layer and inviscid-flow 
equations is obtained. This latter approach is often referred to as interacting boundary- 
layer theory. In cases in which massive separation occurs, such as a stalled airfoil, it 
seems apparent that the solution of the Navier-Stokes equations must be used to obtain a 
detailed description of the flow field. In this case, the separated viscous region occupies 
a significant portion of the flow and thereby makes the inherent assumptions of interacting 
boundary-layer theory highly questionable. Alternately, there are many flows of interest 
in which, despite the presence of separation, the viscous flow is confined to a relatively 
thin layer adjacent to the body surface, thereby rendering the use of the Navier-Stokes 
equations unnecessary for most, if not all, of the flow field. In the present paper, a 



technique is presented for obtaining approximate solutions for boundary layers in which 
the separated region is of a limited extent. Despite this simplifying feature, conven- 
tional or direct boundary -layer theory (pressure distribution prescribed) cannot be used 
without significant modification, such as that used in reference 1, because the boundary- 
layer solution is singular at separation. 

In a recently developed boundary-layer theory, the displacement thickness is pre- 
scribed in the streamwise direction, thereby resulting in the calculation passing smoothly 
through the separation point. The surface pressure distribution is then deduced from the 
boundary-layer calculation. This idea of an inverse boundary -layer procedure to avoid 
the separation singularity was first presented in reference 2 and later was reformulated 
in reference 3. An alternate inverse boundary -layer procedure is to prescribe the skin- 
friction distribution as demonstrated in references 3 to 5 and others. However, the pre- 
scribed displacement-thickness technique is preferable when the simultaneous interaction 
of the boundary layer and inviscid flow is considered. 

The purpose of the present paper is to give a new prescribed displacement-thickness 
procedure, which is considerably simpler than that given previously in reference 3. The 
new procedure is no more complicated than a direct boundary -layer calculation. This 
procedure is presented for incompressible flow, and comparisons are made with the tur- 
bulent experimental data presented in reference 6. With a few modifications, this method 
is then extended to compressible flow, and some sample calculations for laminar and tur- 
bulent separated flow are presented. 

An axisymmetric version of the incompressible procedure has been successfully 
interacted with an inviscid program described in reference 7 in which the relaxation solu- 
tion of the potential equation is obtained. Calculations for the low-speed separated flow 
over an axisymmetric boattail-sting juncture are discussed in reference 8. 

SYMBOLS 




coefficients in finite-difference equations 




coefficients in recurrence equations 


Cj skin-friction coefficient 




skin-friction coefficient multiplied by 



2 



f 

f 

F 

h 


transformed stream function 
perturbation stream function in in-method 
nondimensional velocity ratio, u/u^ 



0 for two-dimensional flow 


1 for axisymmetric flow 


k 

K 


constant of proportionality used in similar form of equations for m-method 


ratio of two adjacent rj-mesh increments 


= PM.. 
Pe^e 


reference length 


m 


perturbation mass flow 


M 


Mach number 


N 


maximum number of points across boundary layer 


N 


Pr 


N 


Prandtl number 


turbulent Prandtl number 


Pr,turb 

q column iteration counter 


Tg body radius 


R 


free-stream Reynolds number, 


p U L 

^oo oo 

M«, 


3 


T 

u 

U, 

V 

V 

X 

y 

a 

ii 


nondimensional static temperature, T ^ 
velocity component parallel to surface divided by 
free-stream velocity 

velocity component normal to surface divided by and multiplied by 

transformed normal -velocity component 

coordinate along surface divided by L 

coordinate along surface divided by L and multiplied by 


= (y - 1)M^ 


pressure-gradient parameter for Levy-Lees formulation, 


pressure -gradient parameter for m-method. 


d(ln u^) 




pressure-gradient parameter for 5 -method, 6 


.2.. 


dx 


Ue d^ 


y 

6 


ratio of specific-heat coefficient 

transformed displacement thickness multiplied by 


displacement thickness 


Jo V PeU^el 


dy 


dimensional displacement thickness 


6F,6u,6/3, 

6 ^, 60 , 


increments in dependent variable between two successive column iterations 


Ax grid spacing in x-coordinate 

Arj grid spacing in ?] -coordinate 


4 



At? 

A^ 


grid spacing in t? - coordinate 
grid spacing in ^-coordinate 

eddy-viscosity coefficient 




N-n 

1+1 


M N 


Pr,turb 


transformed normal coordinate for 6*-method 


transformed normal coordinate for m -method 






O) 


transformed normal coordinate for Levy-Lees formulation 

static -temperature ratio, 

molecular -viscosity coefficient 

transformed tangential coordinate 

density divided by 

nondimensional stream function 

perturbation stream function in 6*-method 

vorticity 


Subscripts: 


edge of boundary layer 


indices for and ry -coordinates, respectively 


w 


wall 


oo free stream 


ANALYSIS 

Incompressible Boundary-Layer Formulation 

The boundary-layer forms of the nondimensional continuity and x-momentum equa- 
tions for incompressible flow are 

^ + ^ = 0 ( 1 ) 
9x ay 


and 


au au 

U + V — = u 

3x 8y ® dx 


a - au 


9y\ 9y 


( 2 ) 


The normal coordinate y and velocity v have 
For laminar flows, the quantity e, which 
unity, whereas for turbulent flow 


been scaled in the usual manner by 
appears in equation (2), is set equal to 


e = 1 + — (3) 

where e/p is the ratio of the eddy -viscosity coefficient to the molecular-viscosity 
coefficient. The eddy-viscosity coefficient is used to relate the Reynolds stress to the 
velocity gradient in the usual manner; thus 


-pu'v' 


9y 


( 4 ) 


In the present calculations, a two-layer eddy -viscosity formulation has been used and is 
discussed by Cebeci and Smith in reference 9. In the inner region the formulation is 
based on Prandtl's mixing length model along with the van Driest damping factor. In the 
outer region Clauser's velocity defect model is used along with Cebeci's low Reynolds 
number correction. This two-layer eddy -viscosity model has been successfully used for 
many attached-flow calculations with mild pressure gradients; for separated flow it 
serves only as an approximate model. The purpose of the present paper is to develop a 


6 



computational technique for separated turbulent flows; at a later time this method, when 
used to make computations for comparison with separated turbulent data, could conveni- 
ently be used to test new and improved turbulence models. 

A key feature of the present inverse boundary-layer method as well as that given 
previously (ref. 3) is the elimination of the continuity equation by introducing the stream 
function 


u = ^ 

ay 


(5) 


V= 

9x 


( 6 ) 


The combination of equation (5) and the definition of displacement thickness 6* gives the 
asymptotic value of the stream function at the boundary-layer edge 

^ - \ {y - (>*) as y - « (7) 

As a result, a perturbation stream function xp can be introduced as 





(^ - u(y - 6*)) 


( 8 ) 


such that \p ^ 0 as y - °o. In this manner a boundary condition is imposed on the 
stream function at the outer edge provided that the displacement thickness is known. 

This condition replaces the usual outer boundary condition u — u^ as y — <». If the 
perturbation stream function is introduced into the governing equations (1) and (2) and the 
normal coordinate is conveniently scaled by 5* to give rj = y/6*, then the governing 
equations become 





( 9 ) 


^ - 6 * ^ 


ax 


ax 


\f^4/ + (77 - l)u6’' 


9^ 97? \ 9t7 


(10) 
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where 


^ = 6*^u 


du 
e 

e dx 


( 11 ) 


These equations can be solved either in the usual direct mode with Ug given or in the 
inverse mode with 6* prescribed and u^ and unknown. In the direct mode 6*, 
which appears in equations (9) to (11), is not the displacement thickness but is some given 
function of x, such as 6* = \[^. In the direct and inverse modes the boundary conditions 
are 

u = ;//=0 at 77 = 0 

Direct; u - u^ as 77 - °c \ ( 12 ) 

Inverse: ip ~ 0 as 77 — 00 

J 

In the inverse mode it is assumed that the y-derivatives vanish at the outer boundary; 
consequently, the x-momentum equation becomes 

6 * 2 u ^ = I ( 13 ) 

ax 

The discussion on the numerical scheme shows that the unknown pressure gradient (3 
is deduced implicitly along with the solution for u and \j/. Hereafter this formulation, 
combined with the numerical procedure described in the following section, is referred to 
as the 6*-method. 

In the previous 6*-method presented in reference 3 , the unknown pressure gradient 
was eliminated by differentiating the x-momentum equation with respect to y, thus intro- 
ducing the vorticity transport equation. The present formulation is preferable because, 
not only can it be solved in either the direct or inverse mode, but it avoids the complica- 
tion of having to solve for the unknown surface vorticity as was required in the previous 
formulation. Furthermore, the vorticity formulation introduces the second y-derivative 
of the eddy-viscosity coefficient, as shown in reference 10 . This formulation was some- 
times found to cause oscillations in the outer part of the boundary layer when the two- 
layer model described previously was used. 
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Numerical Procedure 


The governing equations are solved by using the Crank-Nicolson finite-difference 
scheme with Newton linearization. In the reversed flow region the Reyhner and Fliigge- 
Lotz (ref. 11) approximation 

u — = 0 for u<0 (14) 

is used to prevent instability while preserving the usual rapid forward-marching scheme 
used for the boundary -layer equations. In a previous paper, Carter and Wornom (ref. 12) 
showed that a similar approximation 

u — = 0 for u<0 (15) 

dx ^ ’ 

introduces negligible errors in the reversed flow region if the absolute value of the 
reversed flow velocity is less than approximately O.lOu^. Cebeci (ref. 13) repeated some 
of the calculations made by Carter and Wornom (ref. 12) and showed that the approxima- 
tions given in equations (14) and (15) yield about the same result. 

Finite-difference equations .- A variable grid is used across the boundary layer to 
account adequately for the two-layer structure of a turbulent velocity profile with as few 
grid points as possible. With n used to denote the grid-point index in the normal direc- 
tion, the first and second derivatives are approximated by 


H = ^n+1 ~ %-l 

(1+K)At)j^_^ 


(16) 


and 





n 


2 

K(1 +K)A772_j 



n+1 








(17) 


where “ ’In " ^n-1 the constant K = which is a variable grid that 

has been used by a number of investigators of turbulent boundary layers. The finite- 
difference forms of the governing boundary -layer equations were obtained by using the 
difference expressions given in equations (16) and (17) in conjunction with the Crank- 
Nicolson scheme. Newton linearization was used on the nonlinear terms to accelerate the 
convergence of the iterative solution made at each streamwise location. The linearized 
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finite -difference form of the stream function and x-momentum equations, given in equa- 
tions (9) and (10), are written as 




(6u 
n\ V. 



(18) 


and 




B 6u 
n n 




D 6\p =E 
n n 


Hn5/3 


(19) 


where 6u, 64 /, and 5^ (inverse mode) denote the change in the dependent variable at a 

given point between two successive column iterations; that is 


5u = u^'*’^ - u^ 
n n n 


(20) 


where q is the iteration index. The coefficients in equations (18) and (19), which are 
assumed known from the previous iteration, are given in appendix A. 

The truncation error in equations (16) and (17) is o|(l - K) thus, the scheme 

is second-order accurate provided that K = 1 + 0(A?7). Blottner (ref. 14) used the same 
finite -difference expressions and demonstrated numerically that they are second-order 
accurate. Blottner reports to have made these calculations for a fixed value of K = 1.82; 
however, it can be shown that this statement is in error, and the actual K used by 
Blottner depends on the number of points N used across the boundary layer and is given 
by 


K = 1.82‘»/N 


( 21 ) 


Equation (21) reduces to 

K = 1 + In 1.82 (22) 

N 

as N - “c and thus shows that K = 1 + 0(Ar?), thereby resulting in a second-order 
accurate scheme. 

Newton linearization of viscous term.- The number of column iterations in the tur- 
bulent calculations could be halved by using Newton linearization on the viscous term in 
the inner region. When the shear is positive in that region (the following procedure was 
used only in that case), the viscous term is written as 
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au\ _9_ 

9??\ 9ry 

where Cj is assumed known. With the difference expression given in equation (17) and 
Newton linearization, the viscous term can be expressed as 


1 +c, 

^ ErqJ dri 


( 23 ) 


dv\ d7]) 


2 

K(1 +K)Ar?2_^ 




6u, 


n+1 




X + K/2e 


n-i 
2 


5%_i + e 


n + - 
2 


1 Vl 


(e 1 + Ke 

n--j 


+ Ke 

n 


1 

2 



(24) 


Comparison of equation (24) with the usual formulation, where Newton linearization is not 
used on the eddy -viscosity coefficient 


377 



2 

K(1 + K)A?72_^ 



6u 


n+1 






+ K? i6u^_j + e e ^ + Ke ^ u^ 

2 2 \ 2 2/ 


* n 

2 


(25) 


shows that Newton linearization of the viscous term adds only a slight complication. 

Keller and Cebeci (ref. 15) were the first to use Newton linearization on the boundary- 
layer equations; however, they did not use it on the turbulent viscous term. Nevertheless 
they claimed that their calculations converged quadratically; that is, the residuals decrease 
as 


11 


2 


(26) 


5u 


q+1 


6u. 


with each subsequent iteration. 

This convergence rate was observed only in the laminar calculations in the present 
study; even with the Newton linearization on the inner eddy -viscosity law, the present 
turbulent calculations did not converge quadratically. On the basis of the observations 
in the present paper, it is not clear how Keller and Cebeci (ref. 15) obtained quadratic 
convergence for their turbulent calculations. 

Solution of linearized difference equations. - Repeated application of equations (18) 
and (19) across the boundary layer results in a system of block tridiagonal linear equa- 
tions, which are solved by the recurrence relations 


and 


6^n = Qk - 


(27) 


(28) 


where 






A 






n q^ 


q, = B„ - + D„pJa' 1 - D„(r' , + p) 

^1 n \ n n n/ n-1 n\ n-1 n/ 


(29) 
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and 




Rn - -A;Q2 


^^2 " *^n-l + ^n(^ + ^n-l) 


(30) 


The recurrence constants given in equations (29) and (30) are computed from the wall to 
the edge of the boundary layer. Since the surface boundary conditions are all homo- 
geneous 6uj = = 0^, the recurrence constants are initialized at n = 1 by 

setting 

E\ = H’l = A’ = Q\ = s; = r; = 0 (31) 

In the inverse case the perturbation in the edge velocity 6Uj^ (where N denotes the 
grid point at the boundary-layer edge) and the perturbation in the pressure gradient 6^ 
are deduced at the edge of the boundary layer. These quantities are found from the 
simultaneous solution of the recurrence relations, which results in 

= ^N-1 - (32) 

^’^N-1 (33) 


and the finite -difference forms of equation (13) and equation (9), where the latter is eval- 
uated at n = N - 1/2, gives 

(34) 


and 


13 



in which the boundary condition 5\p^ = 0 has been imposed. 
The solution of equations (32) to (35) results in 

Qn + Qn- 1 ^N^N-1 ^N^N-l) 

5Uj^ = — ^ 

^N-1 + ^N-l) ®n(®N- 1 ^N^N-l) 


and then 6Uj^ j, 54 /-^ p and 6^ are deduced from equations (32) to (34). Now the 
back- substitution process is initiated, and the process continues down to the wall by using 
equations (27) and (28). Note that in this procedure the pressure gradient ^ is deduced 
simultaneously along with the boundary -layer solution u and This procedure is 
similar to that developed previously by Reyhner (ref. 16). In a direct calculation, since 
Ug and ^ are known, 6Uj^ = 6^ = 0 , and the back-substitution process begins by using 
equations (32) and (33). In either the direct or inverse calculation, the column iterative 
procedure is continued until 


where is set equal to 10“'^ in the present calculations, 

tions typically required is 5. 

Compressible Boundary-Layer Formulation 

The compressible formulation is obtained by using a suitable compressibility trans- 
formation followed by an additional transformation to obtain an inverse formulation. 

These transformations could be considered together, but it is simpler to treat them sepa- 
rately. In the present study the Levy-Lees compressibility transformation (ref. 17), is 
used, although an analogous inverse formulation would have resulted if the Stewartson 
(ref. 18) transformation had been used. The nondimensional boundary-layer equations 
expressed after the Levy-Lees transformation (ref. 19) are given as follows; 

Continuity; 



(37) 

The number of column itera- 


2^^+F+^ = 0 

9^ a)) 


(38) 
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Momentum: 


2^F ^ + V — = 

dij 




Energy: 


2|F — + V 

3^ 


9^7 



+ 


^ 

^Pr 9^ 



(39) 


(40) 


where the coordinates are transformed by 

i = j: Pe“e''e4' ** 

and 


p^u^ri, r>y 
77 _ e e B 1 

\[^ ^0 P 




and the dependent variables are given by 


(41) 


(42) 




Additional quantities in equations (38) to (40) are given by 


= PP 
Pe^e 


a = (y- 1)M^ 


/3 




(44) 
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The boundary conditions are given as follows: 


Wall boundary: 

F(^,0) = V(^,0) = 0 ^ 

0 (^, 0 ) = 0 ^( 1 ) 

) 


90 



M (^) 




(45) 


Edge conditions: 

F and 0 — 1 as 77 — <» (46) 

The eddy -viscosity formulation used is the same as that in the 6*-method. The tur- 
bulent Prandtl number enters the analysis in the compressible case in the term 


e 


1 + 


e 

M N 


N 


Pr 


Pr,turb 


(47) 


where 


^ u’v’(9T/9y) 
Pr,turb vtx-(au/9y) 


(48) 


Harris (ref. 19) presents a discussion of the turbulent Prandtl number; in the present 
analysis Np^ = 0.95. 

The development of the inverse compressible formulation proceeds in an analogous 
manner to that developed previously for incompressible flow, A stream function f is 
introduced as 


and the continuity equation is integrated to give 
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( 50 ) 


where the boundary condition V(|,0) = 0 has been imposed. Introduction of a perturba- 
tion stream function is desired so that an outer boundary condition can be imposed from 
the prescribed displacement thickness. This boundary condition replaces the usual con- 
dition of (3, the pressure gradient, being prescribed. The edge value of the stream func- 
tion is 



The compressible displacement thickness is given by 



which becomes, after the Levy-Lees transformation 


.JE-C 


(0 - F) d?7 


Combining equations (51) and (53) results in 


pOO 

= \ 0 d77 

>’o 


Pe^e^’B^* 




A perturbation stream function is defined as 


I 


f = f - F( \ e dr] - 
0 






such that the stream -function boundary conditions are 


f(^,0) = f(|,°o) = 0 


(51) 


(52) 


(53) 


(54) 


(55) 


( 56 ) 


It is convenient to prescribe the perturbation mass flow 
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(57) 


instead of the displacement thickness because this grouping of quantities appears in the 
definition of the perturbation stream function. An inverse boundary -layer procedure 
based on mass flow instead of displacement thickness should be useful in a complete 
viscous-inviscid interaction since the effect of the boundary layer on the inviscid flow may 
be represented by imposing an injection (or suction) velocity at the solid surface given by 


V - ^ 

^ i dx 

Pe^B 

It is also convenient to simplify the normal coordinate by introducing 


(58) 


77 = 

m 

which in terms of physical variables gives 



(59) 


Note that rj is the usual Dorodnitsyn variable (ref. 20) scaled by the displacement thick- 
ness. Use of this variable helps to reduce the effect of boundary -layer growth in the 
computation. With these definitions, the perturbation stream function becomes 

f - f - ^ (n - 1 + h) (60) 

V2| 


where 


h(|) s \ {e - 1) dfj 
-'0 

The inverse compressible formulation is deduced from equations (38) to (40) by 
transforming from t) to t) and from V to f. The resulting system of equations is 
referred to as the prescribed perturbation mass-flow formulation and is given as 
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( 61 ) 


87 ] 




8F 

87 ] 


-2„ 8F - 8 
m F m — 


9^ 


94 


+ mF(7i - 1 + h) 


8F 

dij 



and 


(62) 


— 88 
m F — 
94 



D 21 f + mF(8 - 1 + h)| 


97 ] 



(63) 


- d In uj 

These equations are solved for F, 0, f, and /3 = ^ ^ 

distribution m subject to the following boundary conditions: 


for a prescribed streamwise 


Wall boundary: 

F(4,0) = f(|,0) = 0 
0(4,0) = 8^(4) 

or I 



Edge conditions: 


(64) 


F = 8 — 1 and f — 0 as t] — <» (65) 

Hereafter this formulation, combined with the numerical technique, which is now described, 
is referred to as the m-method. 

The Crank-Nicolson finite -difference scheme with Newton iteration is used to solve 
equations (61) to (63). The resulting finite -difference equations and solution technique are 
presented in appendix B since the technique is the same as that described previously. The 
energy equation is uncoupled from the combined stream function and x-momentum equa- 
tions; as a result, quadratic convergence is not observed in the iterative convergence. At 
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a later time, when the present scheme is embodied in a complete viscous-inviscid inter- 
action problem, it is planned to treat all three governing equations in a coupled manner 
by using Newton iteration. The resulting finite-difference equations will be a block tri- 
diagonal system with each block being a 3 x 3 matrix provided that h(|) is assumed to 
be known from the previous boundary -layer sweep. 

RESULTS AND DISCUSSION 
Comparison With Experimental Data 

Chu and Young (ref. 6) reported an experiment in which measurements were made 
on a flat plate approximately 2^ m long, which was first subjected to a favorable pressure 
gradient and then to an adverse pressure gradient. This pressure distribution was 
obtained by placing a 7. 5-cm -diameter circular cylinder 11 cm above the flat plate with 
its axis parallel to the plate and normal to the windstream. 

Calculations were made for this low-speed flow (u^^ = 18.7 m/s) with the 6* 
deduced by Chu and Young from velocity profile measurements prescribed in the stream - 
wise direction. The 6*-method is used instead of the incompressible version of the 
m-method since, with the latter formulation, the computation is done in terms of a trans- 
formed x-coordinate 

^ oo 

^ Ug dx (66) 

and requires m = u^6* to be prescribed as a function of Since u^ is unknown, it is 
preferable in the present application of inverse boundary -layer theory to use the 6*-method 
which requires only 6*(x) to be specified. Note that in a complete interaction calculation, 
in which an appropriate inviscid scheme is combined with the m-method, repeated stream- 
wise computations are made, and therefore the most recent distribution of u^ can be 
used to transform from the x-coordinate to the ^-coordinate and conversely. 

The upstream profile of u was obtained by using Coles' (ref. 21) wall-wake velocity 
distribution. This profile formulation contains three constants, which are determined by 
matching with experimental data. In the present case four quantities, 6*, C|, 6, and Ug, 
were known from Chu and Young's measurements. In order to use all the data to deduce 
the constants in Coles' formulation, a least-squares error calculation was used. Also, 
the wall-wake profile was modified near the surface by Kleinstein's formulation (ref. 22), 
thereby accounting for the laminar sublayer. At the outer edge of the boundary layer an 
exponential function was used to force the wall-wake profile for u to approach asymp- 
totically Ug. After the u-profile was known, the perturbation stream function {p was 
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deduced by integrating equation (9) from rj = 0 where = 0. Ninety-three grid 

points were used across the boundary layer with K = 1,09. 

A cubic spline fit to the displacement-thickness distribution measured by Chu and 
Young is shown in figure 1. This curve fit to the data was obtained by using a computer 
program developed by Smith and others (ref. 23) to represent the desired function by a 
cubic polynomial over prescribed intervals. Note that the displacement thickness shown 
in figure 1 is dimensional and is denoted as 5*. In order to demonstrate the capabilities 
of the present technique in separated flow, the prescribed 6*-distribution extends in fig- 
ure 1 beyond the data of Chu and Young in a manner typical of a turbulent boundary layer 
incurring both separation and reattachment. 

Figure 2 shows the comparison of the computed skin-friction coefficient on the flat 
plate with the experimental distribution obtained by Chu and Young. The agreement is 
good despite the use of an algebraic eddy -viscosity model, which is generally considered 
to be inappropriate for turbulent boundary layers subjected to large pressure gradients. 
Chu and Young also applied several direct boundary -layer procedures to these data and 
found that only the Bradshaw (ref. 24) and Kuhn and Nielsen (ref. 25) procedures gave 
good agreement with the data. None of these direct procedures can be used downstream 
of the separation point as exhibited by the present inverse calculation shown in figure 2. 
Figure 3 shows good agreement between the computed and measured velocity distribution 
at the edge of the boundary layer on the flat plate. Although more data are needed to 
assess the present method in a separated region, the good agreement shown in figures 2 
and 3 between theory and experiment is encouraging and motivates the further use of 
inverse boundary-layer theory. The present comparison was simplified by the use of 
the experimental 5* -distribution, thereby requiring only one streamwise boundary -layer 
calculation to be made. In a complete interaction calculation it is necessary to combine 
a suitable inviscid calculation with an inverse boundary-layer calculation and to solve the 
two simultaneously until convergence. A study demonstrating this technique for incom- 
pressible flow is reported in reference 8. 


Results Obtained With m-Method 

Similar solutions.- A similar form of the mass-flow prescribed formulation is 
obtained by setting m = k^Y^, where k is a prescribed constant, and by assuming that 

— =0. The similar form of these equations for incompressible laminar flow is given by 


r = k(i - ?7) F' 


(67) 


and 
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-kfi + k(7j - 1 )f]f’ = k^;3(l - F^) + F" 


( 68 ) 


with the boundary conditions 
f(0) = F(0) = 0 


(69) 


and 


f — 0 and F ^ 1 as 77 — 00 (70) 

where the prime denotes differentiation with respect to fj. A general, compressible, 
similar formulation cannot be obtained from the m-method because the edge Mach num- 
ber Mg appears explicitly in the energy equation. H /3 / 0, then Mg = Mg(|), and thus 
the energy equation is not independent of A compressible similar formulation can be 
obtained by using the total -enthalpy form of the energy equation as was done by Cohen 
and Reshotko (ref. 26), thereby eliminating Mg from the equation when the Prandtl num- 
ber is unity. 

Incompressible similar solutions were obtained by solvii^ equations (67) and (68) 
for different values of k subject to the boundary conditions given in equations (69) 
and (70). The prescribed constant k is easily shown to be equal to 6, where 6 is the 
transformed displacement thickness expressed in terms of the Levy-Lees coordinate 77 ; 
thus 

6 = (1 - F) d77 (71) 


Results of these solutions are compared in figure 4 with so-called exact solutions, which 
previous investigators (e.g., see ref. 27) have obtained by numerical means. Excellent 
agreement is obtained. It is observed in figure 4 that for each value of /3, the pressure- 
gradient parameter, there are two values for both the transformed displacement thick- 
ness 6 and the wall shear F^. Thus, in a direct calculation, that is when /3 is 
prescribed, it is not clear what values of F^ and 6 will be obtained. Previously, 
unless a shooting technique was used, the forward flow solution F^ > 0 and corre- 
sponding value of 6 always resulted for a given value of p. Werle and Bertke (ref. 28) 
presented a double iterative procedure to obtain the reverse flow solutions F^ < 0, which 
were originally discovered by Stewartson (ref. 29). In the present approach there was no 
difficulty in obtaining either the forward or reversed profile because figure 4(a) indicates 
that, given a value of 6, a single value of /3 results. 
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Laminar separation.- As a demonstration of the present procedure for computing 
separated boundary layers, the m-distributions shown in figure 5(a) as a function of ^ 
were prescribed for both an incompressible and a compressible calculation. In the com- 
pressible case, = 2 and the wall temperature was set equal to the free -stream 
stagnation temperature with Np^ = 1. The two m-distributions in figure 5(a) differ by a 
constant so that both m-distributions match that of a flat-plate solution at ^ = x = 1. 

The initial boundary-layer profiles at this point were flat-plate similar solutions. 

The resulting calculations for both cases are shown in figures 5(b), 5(c), and 5(d) 
for the edge velocity, skin friction, and displacement thickness, respectively. These 
quantities are plotted against the physical coordinate x, which is obtained after the cal- 
culation is completed from 


Jt^n P„u, 




which for isentropic flow and a linear temperature-viscosity law becomes 




AL 


1 + ^ ^ (] 

2\ 

^ J 

2 “\ 



y/(y-l) 


(72) 


(73) 


Although figures 5(b) to 5(d) show that the calculations for = 0 and = 2 differ 
in detail, the resulting solutions show the same basic features. This result is expected 
because equation (73) shows that if u^ remains close to unity (as shown in fig. 5(b)), 
then the physical and transformed streamwise coordinates are nearly the same, independ- 
ent of the free-stream Mach number. Both solutions are smooth, and no difficulty was 
encountered near the separation or reattachment point in either solution. Asa check on 
the computer program for the ra-method, the displacement thickness shown in figure 5(d) 
for the incompressible case was prescribed for the 6*-method. The resulting solutions 
for Ug and c^ were identical to those shown in figures 5(b) and 5(c), respectively, 
for M =0. 

oc 

Turbulent separation. - A turbulent computation, similar to that shown for laminar 
flow, was made to demonstrate the present procedure for both an incompressible and 
compressible separated turbulent boundary layer. The prescribed m-distributions are 
shown in figure 6(a) for the calculations with = 0 and = 2. As was the case for 
the laminar calculation, the prescribed m-distributions match a turbulent flat-plate solu- 
tion at the initial streamwise station x = ^ = 0.95 ^ = 0.95 x 10®^ When this cal- 

culation was made the computer program using the m-method had not been written so that 
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a direct calculation could be made. Hence, the initial boundary- layer profiles had to be 
obtained elsewhere. In the incompressible case the program for the 6*-method was run 
in the direct mode with 3 = 0, and the flow was assumed to be turbulent from the leading 
edge of the flat plate. For the compressible case the upstream profiles were obtained 
with the computer program reported in reference 30. 

The resulting distributions of edge velocity, skin friction, and displacement thick- 
ness are plotted against the physical coordinate x in figures 6(b), 6(c), and 6(d), respec- 
tively. In both cases, flow separation and reattachment occurred, and no difficulty was 
encountered in obtaining a smooth solution. In contrast with the laminar case, the distri- 
butions of Ug shown in figure 6(b) differ significantly from unity, thereby resulting in 
different inverse transformations from the |-coordinate to the x-coordinate using 
equation (73). 

In view of the different distributions of Ug, it is not surprising that large differences 
occur between the distributions of skin friction and displacement thickness for = 0 
and = 2 shown in figures 6(c) and 6(d), respectively. A common feature is that both 
solutions show that the points of reattachment and maximum displacement thickness occur 
close to the point of minimum-edge velocity (maximum pressure). This trend has been 
observed in separated turbulent flow measurements, such as those of Alber and others in 
reference 31. 

Viscous interaction . - With the development of the m-method described in this paper, 
efforts should now be directed at the complete problem of viscous interaction. This 
inverse boundary -layer procedure, which to the author's knowledge is the only compres- 
sible, inverse, finite -difference procedure based essentially on displacement thickness, 
appears to be well suited for viscous interaction as shown by equation (58). As mentioned 
previously, a study has been completed using the 6*-method on axisymmetric boattails and 
is reported in reference 8. The results and interaction scheme seem promising. Several 
compressible interaction calculations using inverse integral procedures have been 
reported in the literature, such as those by Kuhn and Nielsen (ref. 4) and Thiede (ref. 32). 
Integral techniques have the disadvantage that the shape of the velocity profile must be 
specified by the user, whereas the finite -difference calculation does not. Also, the present 
finite -difference inverse boundary-layer procedure should provide a good means to test 
various turbulence models for separated flow much more economically than can be done 
with solutions of the Navier-Stokes equations. Most of the current Navier-Stokes compu- 
tations are forced to use relatively coarse grid spacing in the viscous region adjacent to 
the body region. In constrast, boundary -layer procedures provide complete resolution of 
the viscous region because the computation is confined entirely to the viscous region. 

The present m-method should be useful for separated flows with supersonic as well 
as subsonic edge conditions. The main motivation for the development of inverse 
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boundary-layer theory has been directed at subsonic separated flows to eliminate the 
singularity encountered at the separation point. Werle and Vatsa (ref. 33) developed a 
new supersonic viscous interaction technique, which was a major achievement because 
it converted the calculation from an ill-posed initial value problem to a well-posed 
boundary-value problem. The disadvantage of their procedure was that, unless a great 
deal of care was taken, a singularity was obtained at the point of separation. If the pres- 
ent ih-method were used, while retaining the interaction model of Werle and Vatsa, it is 
anticipated that a well-posed boundary-value problem would result with no singularity at 
the separation point. 


CONCLUDING REMARKS 

Two inverse boundary -layer formulations and solution techniques are presented for 
separated boundary -layer calculations. In the first technique, which is only for incom- 
pressible flow, the displacement thickness is prescribed; comparisons made with this 
method show good agreement with experimental low-speed turbulent data up to the point 
of separation. Further comparison with experimental data is needed to assess the 
accuracy of the present boundary -layer method in a region of separated flow. It is antici- 
pated that the boundary -layer approximations are still valid provided that the separated 
region is not large but, based on the current literature, significant improvements are 
needed in the turbulence model for a separated region. 

The second inverse technique presented in this paper is applicable to compressible 
flows and is based on a prescribed perturbation mass-flow distribution. A solution tech- 
nique is presented for this formulation, and sample calculations are presented for laminar 
and turbulent incompressible and compressible flows. In all these calculations the solu- 
tion passes smoothly through the points of separation and reattachment. 

Both inverse procedures appear to be well suited for use in complete viscous- 
inviscid interactions because they are based primarily on the displacement thickness. 

The next step in the development of a complete viscous-inviscid interaction theory is to 
combine the compressible inverse procedure with an appropriately modified inviscid tech- 
nique for subsonic or supersonic flows. This viscous-inviscid iteration technique can 
then be compared with experiment to test its accuracy as well as the simultaneous 
development of more sophisticated turbulence models. Interacting boundary-layer theory 
provides a more economical means by which different turbulence models for separated 
flows can be tested than current Navier-Stokes computations. Finite-difference procedures 
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should prove superior to integral techniques because the latter requires the specification 
of the shape of the velocity profile. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
July 19, 1978 
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APPENDIX A 


COEFFICIENTS OF FINITE -DIFFERENCE EQUATIONS FOR 6*-METHOD 


The stream function and x-momentum equations, given by equations (9) and (10), are 
differenced about the points ^m, n - and respectively, where m denotes 

the x-grid-point index and n denotes the r? -grid-point index. The resulting finite - 
difference equations are given by equations (18) and (19) in which the following coefficients 
appear: 


An = e - 


m,n - - 
' 2 

(1 +K)A7?2_^ 


(Al) 


B„ = 2au - fc, + 
n m,n ’ 


^ IH-Kc 1 

m,n+- m,n-- 

2 ’ 2 


^ K(1 + K) Avl_^ 


(A2) 


Cn - -e - 


m,n + - 
2 


K(1 +K)At7^_j 


(A3) 


Dn = -bif 


(A4) 


^ 1 + a[u^ 1 - u^ I + ef 

" m-i \ m-l,n ^m,n/ + 


K(1 +K)At72_^ 


„^l“m,n+l"[^ 1+Ke i)Vn 

,n + - ’ V ""’"+2 ‘"’"" I / ’ 


c 1 u 
m 


+ Ke 


m 


" 'm-l,n.i"-l>--l " rm-l,n.i " Vl,n 

2 2 / 


+ Ke 1 u 1 , 

m-l,n-i m-l,n-l 
’ 2 


(A5) 
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2 Ax 


(A6) 

(A7) 

(A8) 

(A9) 



(AlO) 


(All) 


(A12) 


(A13) 


e = bi!i' -bot^ 1 +CiU - Cf,u , 

l^m,n 2^m-l,n 1 m,n 2 m-l,n 

and 

^ “ ^m,n+l ” ^m,n-l ^m-l,n+l ^m-l,n-l 


(A14) 


(A15) 
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Note that 



is given by 




(A16) 
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APPENDIX B 


FINITE-DIFFERENCE EQUATIONS AND SOLUTION TECHNIQUE 

FOR m -METHOD 


Application of the Crank-Nicolson scheme with a variable grid and with Newton 
iteration to equations ( 61 ) and ( 62 ) about the points ^m, n - and . respec- 

tivelyj results in the following system of linearized finite— difference equations where the 
notation is the same as that used previously: 


(Bl) 


«n-“„-l + Pn(^fn-'>Vl)=«n 


where 






m,n-- 


(1 + K)AT)2_^ 


Bn - 2aF^,n ' ^1^2 j 


K(1 + K) Arj2_^ 


(fe) 1 + K(fe) j 

ni,n + - 


J 


m 


^n=-A - 


m,n + - 


K(1 + K) A7?2_^ 


(B2) 


(B3) 


(B4) 


(B5) 


Bn " “^1^2 


(B6) 


30 



APPENDIX B 


“ ^m,nj 1 

' m -- 


ile) 


K(1 +K)Af;2_^^"''-^n,,n + i 


Ue) + K(Ze) . 

m,n+- m,n-i 


+ K(Ze) 1 n 1 + 1 F 1 1 

m,n-i m-l,n + l m-l,n+l 


ile) j + K(Zi) 


m-l,n 


+ K(Ze) . F 1 ^ 

m-l,n - - ni-l,n-l 


Hn=m" 1 
*"-2 


^m,n ^ ^m-l,n _ | ^m,n ^ 


Pn = 





Q 


n ^m,n-l ■^m,n + ^n(^m,n-l “ ^m,n) 


The additional coefficients which appear in equations (B3) to (BIO) are given by 


a = 


m2 

m -- 
2 

2A| 


F 

m,n 


(B7) 


(B8) 


(B9) 


(BIO) 


(Bll) 
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m 


b, = 




1 2(1 + 2^J j 


m - ■ 


(B12) 


m 


bo = 




2 2(1 + K) \ A| 1 


m-- 


(B13) 


1 




(B14) 




(B15) 


®1 " ^l^m,n " ^2^m-l,n ^l^m,n ” ‘^2^m-l,n 


(B16) 


^2 ^m,n+l 


^m,n-l ''' ^m-l,n+l 


- F. 


m 


-l,n-l 


(B17) 


In the present calculations for simplicity I = = 1, which is an approximation fre- 

be^e 

quently used by others. Also, Newton linearization of the eddy viscosity in the 
x-momentum equation was used in the same manner as described in the main text for the 
6*-method. 


As in the 6*-method equations (Bl) and (B2) can be solved either in the inverse mode 
where m is prescribed and ^ is deduced or in the direct mode where ^ is prescribed. 
In the direct mode if m = then 77 = !), = / 3 , and the resulting formulation is the 

usual Levy-Lees boundary -layer formulation where the transformed normal-velocity com- 
ponent has been replaced by 


r 1 


v= -72li<!/2l] 

|1 + F(/7 - 1 + h^j 


(B18) 
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Equations (Bl) and (B2) are solved by the same solution technique presented for the 
incompressible formulation in the main text. In the inverse case, 6/3, 6Fj^ j, and 
^^N-1 deduced from the simultaneous solution of 

~ ^N-1 ” 

= Qn-1 - (B20) 


®^N-1 * ^N^^N-1 


= -Q 


N 


(B21) 


where the boundary conditions = 0 have been applied. In the direct case equa- 

tion (B21) is replaced by 


“ '^^N-1 “ ^N^^N-1 " 

and equations (B19), (B20), and (B22) are solved for 6F^ 6f^ and 6f^ with 

6/3=0. The back-substitution procedure is then carried out from the boundary-layer 
edge to the surface. 

After F and f are updated by solving the x-momentum and stream-function 
equation, the linearized form of the energy equation is solved. The finite -difference form 
of this equation is given by 




where 


m 


K = ^i- 


m.n — 
’ 2 


Npj,{l +K)At/^_j 


- ^(^m,n + ^m-l,n) + 


in) 1 + K(Ze) 1 
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Np^K(l +K)Ar/2_j 


(B23) 


(B24) 


(B25) 
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The solution to the energy equation is given by 


“n = - A„6e„.i 


(B28) 


where 


D - C D 1 
- n n n+l 
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The constants and are evaluated from the outer boundary n = N with 

~ ® wall. If the wall is isothermal, then 60 j = 0 and the solution for 

60j^ is deduced from equation (B28). If the heat transfer 90/&0{^,O) is specified, then 
the finite-difference form of this Neumann condition is solved along with equation (B28) 
to deduce 60^. Then the back- substitution process given in equation (B28) begins. 
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Figure 2.- Skin-friction distributions. = 2.08 x 10®; = 18.7 m/s 
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